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Abstract 

We make a review of the two principal models that allows to explain the imbibition of fluid 
in porous media. These models, that belong to the directed percolation depinning (DPD) 
universality class, where introduced simultaneously by the Tang and Leschhorn [Phys. Rev 
A 45, R8309 (1992)] and Buldyrev et al. [Phys. Rev. A 45, R8313 (1992)] and reviewed 
by Braunstein et al. [J. Phys. A 32, 1801 (1999); Phys. Rev. E 59, 4243 (1999)]. Even 
these models have been classified in the same universality class than the Kardar-Parisi-Zhang 
equation [Phys. Rev. Lett. 56, 889, (1986)] with quenched noise (QKPZ), the contributions 
to the growing mechanisms are quite different. The lateral contribution in the DPD models, 
leads to an increasing of the roughness near the criticality while in the QKPZ equation this 
contribution always flattens the roughness. These results suggest that the QKPZ equation does 
not describe properly the DPD models even when the exponents derived from this equation are 
similar to the one obtained from the simulations of these models. This fact is confirmed trough 
the deduced analytical equation for the Tang and Leschhorn model. This equation has the same 
symmetries than the QKPZ one but its coefficients depend on the balance between the driving 
force and the quenched noise. 

1 Introduction 

In the last years there has been a growing interest in the understanding a vast variety of scale 
invariant phenomena occurring in nature. Experiments and observations indeed suggest that 
many physical systems develop correlations with power law behaviour both in space and time. 
However, the fact that certain structures exhibit fractal and complex properties does not tell 
us why this happens. Pattern formation, aggregation phenomena, biological systems, geological 
systems, disordered materials, clustering of matter in the universe and many fields in which 
scale invar iance has been observed as a common and basic feature. 



^E-mail: lbrauns@mdp.edu. ar 
^E-mail: tasio@fcu.um.es 



1 



A crucial point to understand is therefore the origin of the general scale-invariance of natural 
phenomena. This would correspond to the understanding of the origin of fractal structures from 
the knowledge of the microscopic physical processes at the basis of these phenomena. 

In scale invariant phenomena, events and information spread over a wide range of length 
and time scales, so that no matter what is the size of the scale considered one always observes 
surprisingly rich structures. These systems, with very many degrees of freedom, are usually so 
complex that their large scale behaviour cannot be predicted from the microscopic dynamics. 
New types of collective behaviour arise and their understanding represents one of the most 
challenging areas in modern statistical physics. 

The concepts of scaling and power law behavior was introduced to the study of critical 
phenomena in second order transitions. The physics of complex systems, however, is new with 
respect to critical phenomena. The theory of equilibrium statistical physics is strongly based on 
the ergodic hypothesis and scale invariance develops at the critical equilibrium between order 
and disorder. Reaching this equilibrium requires the fine tuning of various parameters. In usual 
critical phenomena the same exponents that define the onset of magnetization also describe 
the liquid vapor transition in water. This strong universaUty appears to be a characteristic of 
equilibrium systems. On the contrary the origin of the scale invariance in nature in the rich 
domain of nonequilibrium systems is not so well understood. Systems far from equilibrium, do 
not seem to exhibit the same degree of universality as the fractal dimension can be easily altered 
by simple changes in the growth process. This lack of universality is sometimes viewed as a 
negative element because one is forced to describe specific systems instead of a single universal 
model. This is why such a problem can only be investigated using many tools as computer 
simulations, analytical tools and suitably designed experiments. While the theoretical activity 
is focused mainly on Monte Carlo simulations, it is very important to understand the relations 
between theory and real experiments. 

Fractal geometry provided the mathematical framework for the extension of these concepts 
to a vast variety of natural phenomena. A principal subject where fractals play an essential 
role is in the study of self-affine features arising in complex systems with many degrees of 
freedom, such as of interface growth in disordered systems including percolation properties of 
fluid displacement in disordered media. 

All these models have been extensively studied by computer simulations. The informa- 
tion generated by these computer experiments brings a visual intuition as a valuable tool in 
science. In this respect computer simulations represent an essential method in the physics of 
these systems, that allow us to design theoretical experiments tested with a computer. Numer- 
ical simulations define the basic characteristics of the models and gives a useful path to their 
theoretical understanding. 

Prom a mathematical point of view the problems explored are particularly difficult because 
they consist of iterative systems with many degrees of freedom and irreversible dynamics. Very 
little can be predicted a priori for systems of this complexity even when sometimes it is not very 
difficult to pose a model that capture the essential features of the phenomena. 

While the great majority of the theoretical activity is based upon toy models it is very 
important to build a bridge between theory and real experiments and this is another basic 
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task of computer simulations. To do that, it is necessary the development of models more 
realist and large scale simulations which can be used also in material characterization. The 
consequence of this approach is the application of fractal concepts to the solution of particular 
experimental problems (Oil industry, disordered materials, phase nucleation, crystal growth 
etc.). The theoretical effort in this field can be separated into phenomenological or scaling 
theories and microscopic theories. At a phenomenological level scaling theory, inspired to usual 
critical phenomena, has been successfully used. This is the nexus to the understanding of the 
results of computer simulations and experiments. This method allows us to identify the relations 
between different properties and exponents and to focus on the essential ones. Moreover, it will 
be useful to gain more insight into the microscopic dynamics that evolves into the macroscopic 
behavior. The connection between the microscopic and macroscopic behavior has not been 
extensively studied. 

In this work we make a review of models of growth in quenched disordered media. In 
Sec. |2| we introduce the dynamics and static scaling used in experiments and models of growth. 
We make also a discussion about the principal source of disorder. In Sec. ^ we make a short 
summary of the main experiments on fluid-fluid displacement. In Sec. § we present the most 
used phenomenological equation that is said to belong to the same universality class than the 
experiments and/or models with quenched media in special to the directed percolation depinning 
(DPD) models. In Sec. ^ we present the Tang and Leschhorn model In Sec. ^ we present 
the Buldyrev et al. model Q. In Sec. some comparisons between both models. In Sec. ^ we 
argue why the Kardar-Parisi-Zhang equation |^ with quenched noise does not describe the DPD 
models even if they are said to belong to the same universality class. In Sec. ^ we deduce an 
analytical stochastic differential equation for one of the DPD models starting from its microscopic 
equation. Finally in Sec. 10 we conclude with a discussion. 



2 Scaling properties of growth in quenched media 

The investigation of rough surfaces and interfaces has attracted much attention, for decades, 
due to its importance in many fields, such as the motion of liquids in porous media, growth of 
bacterial colonies, crystal growth, etc. Much effort has been done in understanding the processes 
that induces the roughness in these fields. When a fluid wet a porous medium a nonequilibrium 
self-afhne rough interface is generated. The interface has been characterized through scaling 
of the interfacial width w = {[hi — (/ij)]^)^''^ with time t and lateral size L. The result is the 
determination of two exponents /? and a called dynamical and roughness exponents respectively. 
The interfacial width follows 

w ~ if t » L"/'^ , (1) 

w^t^ if t < . (2) 

The crossover time between this two regimes is of the order of L"-/^ . 

The disorder affects the motion of the interface and leads to its roughness. The main disorder 
proposed has been the "annealed" noise that only depends on time and the "quenched" disorder 
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due to the inhomogeneity of the media where the moving phase is propagating. While the 
annealed disorder pushes forward the interface the quenched disorder brakes this advancement. 
This last disorder can be due to impurities, fluctuations in the capillary size or it could represent 
all the effects of an inhomogeneous media. 

Some experiments such as the motion of liquids in porous media, where the disorder is 
quenched, are well described by the directed percolation depinning model. Some other exper- 
iments like the flow of fluid trough a disordered media, where the disorder is also quenched 
IQ, |5|, ^, give roughness exponents scattered between 0.6 and 1.25 questioning the existence of 
universality, the foundation of the scaling hypothesis (^. 

3 Experiments on fluid-fluid displacement 

Several fluid-fluid displacement experiments that can be described in terms of simple self-affine 
structures were motivated by the growing interest of physicians on the dynamics of rough sur- 
faces. The experiments are easy to perform but the process are very complex to understand. 
One of the experiments that was first performed in order to explain the fact that the roughness 
exponent a was not predicted by equations with thermal noise was carried out by Rubio et al. 

in 1 -|- 1 dimensions. In this experiment the air is displaced by a dyed water in a porous 
media consisting of glass beads packed randomly into a thin horizontal cell of glass covered with 
Teflon]^. The water was injected by one of its edges wetting the porous media and the fluid- fluid 
interface is digitized. The roughness exponent was measured over a distance i obtaining a value 
of a = 0.73 ± 0.03. This value does not depend on the bead diameter and the capillary number 
Cap0 in the range that they worked. Horvath et al. reanalyzed the interface of Rubio et al. 
and obtained a = 0.81 it 0.08. This discrepancy is not clear, but it could be assign to the sensi- 
bility of this exponent to the details of the data analysis. Also in this experiment the exponent 
(3 ~ 0.65 was obtained. Other value of a = 0.63 it 0.04, where the exponent was measured in the 
saturation, was obtained by Buldyrev et al. [Q, ^ with aqueous fluids absorbed in many kind of 
papers. In these experiments, the scaling value does not depend on the temperature, humidity 
and kind of paper. They inferred that the evaporation is irrelevant. Horvath and Stanley |^] 
performed an experiment in which the evaporation was controlled by confining the paper be- 
tween two transparent polymer sheets. The lower end of the cell was immersed on liquid. The 
shape of the interface was digitized and the mean height of the interface was calculated in real 
time in order to maintain it at a fixed height. This prevents any change in the velocity. They 
found (3 = 0.56 ± 0.03. 

The scattering between the exponents obtained in these experiments puts in doubt the uni- 
versality of their results. Sometimes, the uncertainties could be the result of complex crossover 
behaviors before the asymptotic regime is reached. On the other hand, sometimes, the asymp- 
totic regime is not reached in the experiments. Also in many experiments the interval where 

^DuPont trademark. 

*Cap ~ a^v/r, where a is the diameter of the bead, i; is the velocity of the water and F is the interfacial 
tension. 
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Experiments 


a 


/3 


References 


Flux of fluid 


0.73 
0.81 
0.65 - 0.91 


0.65 




i 


Wetting paper 


0.63 (d=l) 
0.5 (d=2) 
0.62 - 0.78 (d=l) 


0.56 
0.3-0.4 


1 


i 


Bacterial growth 


0.78 




1 




Combustion front 


0.71 









Table 1: Measured exponents from experiments 



the exponents were measured is very short enough (less than two decades), in these cases the 
results obtained are less reliable. This does not mean that the idea of scaling behaviors must 
be rejected, moreover some scaling hypothesis must be reviewed. In table |^ we shows some 
experimental results. 



4 Phenomenological equation of growth in quenched media 

Motivated by the success of the Kardar-Parisi-Zang (KPZ) equation in describing the interface 
motion with thermal noise, it has been proposed that many interfaces in porous media are 
described by the quenched KPZ (QKPZ) equation 

= + vV^h + \{yhf + r?(x, K) (3) 

where r/(x, K) represent the quenched white noise in the media and T is the driving force re- 
sponsible of the advance of the interface. This equation predicts the existence of a depinning 
transition: for driving forces T < Tc the interface is pinned, and for JF > .Fc it moves with 
a velocity v ~ fl^^ where Q is the velocity exponent and /red = [T — Tc)/^c is the reduced 



force. Numerical studies 13 1 indicated that the discrete models can be grouped into two 
universality classes depending on the behavior of the QKPZ nonlinearity A. Isotropic models, i.e. 
models that have no growth direction determined by the random forces, have A = or A ^ as 
fred 0. The scaling exponents can be determined analytically, in one dimension, from Eq. (^) 



with A = 0, obtaining = 1, /3j = 0.75 and 9i = 0.33 [14|. However, anisotropy can induce 
a relevant A at the depinning transition. The exponents characterizing this anisotropic univer- 
sality can be obtained by an exacting mapping to directed percolation, obtaining Oa = 0.63, 
Pa = 0.63 and 9a = 0.63 [Q, ^. Many numerical simulations have been done from a tilted inter- 
face |T2| , ^] in order to justify if the nonlinear term of the QKPZ is present in the dynamics of 



the proposed models. Recently Reka et al. [16| implemented a method that apply to numerical 
models and experimental data and allows to classify the data into one of the two universality 
classes without relying on the determination of the exponents. They assumed that the local 
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velocity has the same scahng behavior than the global velocity v, then u{J-, s) ^ [J- — J-^{s)]^ , 
where u{T, s) is the averaged local velocity of the interface with slope s, ^^(s) is the depinning 
threshold corresponding to these segments, and 6 is the local velocity exponent. In isotropic 
media J^^i^) = Tc and Q = 9, so the ratio V = u{J^, s)/v{T) does not depend on the driven force. 
In contrast, in anisotropic media the depinning threshold, ^^{s) decreases with s and 9=1. 
So, for anisotropic media V has a systematic dependence on Plotting V as a function of s for 
the several models and for the fluid-fluid experiment they concluded that this last experiment 
belongs to the isotropic class. All these numerical results confirms the relevance of a nonlinear 
term for any model but they do not assert to prove if these models are well represented by the 
QKPZ equation or its linear version, the quenched Edward- Wilkinson (QEW) ||T7[| . 

In particular in models with quenched noise belonging to the anisotropic class the macro- 
scopic dynamics is influenced by the coupled effect of the interface itself and the disorder as we 
shall see in Sec. |9|. 

A powerful method to derive the macroscopic behavior of these models is to attempt to 
write the microscopic equation for the evolution of the local height in these models. This is 
done via the master equation approach or directly by writing the microscopic rules for the 



evolution of the local height |2l|]. These two methods have been proved to be equivalent 

m. These microscopic equations that captures the mechanisms of the growth for each model 
must reproduce the macroscopic behavior of the relevant observable that in this case is the 
interface itself. Its continuum version, in a coarsed grained scale is the equation that represent 
its universality class. Let us to introduce the two main models describe the principal features 
of the experiments of fluid imbibition on paper sheet [0]. 



5 Growth mechanisms for the Tang and Leschhorn model 

We present now the microscopic Tang and Leschhorn (TL) model. The interface growth takes 
place in a lattice of N points of edge L (N = L/a) where a is the lattice constant. Usually in 
discrete models a is taken as one, so N = L. We will retain the lattice constant because we will 
use it below. 

Periodic boundary conditions are used. The quenched noise is represented by a random 
pinning force g{r) uniformly distributed in [0, 1] assigned to every cell of the lattice. This 
random force models the random distribution of the fiber of paper. For a given pressure p, 
that represents the driving force that push the fluid into the media, the cells are divided in two 
groups, active (free) cells with g{r) < p and inactive (blocked) cells with g{r) > p. Notice that 
q = 1 — p is the density of blocked cells. It is well known that there exist a critical density of 
blocked cells qc = 0.539, what is the directed percolation threshold for an infinite lattice, above 
qc the interface becomes pinned. 

In this model, the growth event is defined as follow: 

1. If hi is greater than either or /ij+i by one or more units, the height of the lower of 
the two columns (i — 1) and (i + 1) is incremented in a (in case of tie, one of the two is 
chosen with equal probability). 
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2. In the opposite case, hi < min(/ij_i, hi+i) + 2a, the column i advances by one unit provided 
that the cell to be occupied is an active cell. 



3. Otherwise no growth takes place. 

In this model, the time unit is defined as one growth attempt. In numerical simulations at each 
growth attempt the time t is increased by St, where 6t = 1/N. 

Notice that the condition for the interface becomes pinned is that: 

• Ah = zba, and 

• all the sites above the interface are blocked. 

So, the interface becomes pinned if a cluster of inactive sites connected horizontally or diagonally 
expands all the lattice. This happens for q > Qc. This define a cluster of directed percolation. 
For q < Qc the interface is characterized by two correlations lengths that behave approaching to 
the critical point as, 

^11 ~ \q-qcn^ , (4) 
U ~ \q-qcr^, (5) 

which are the characteristic lengths of these clusters in both directions, with z^y = 1.733 and 
u_L = 1.097. 

The interface is specified by a set of integer column heights hi {i = 1, . . . ,N). During the 
growth, a column is selected at random with probability 1/A^ and compared its height with 
those of its neighbors. In a temporal step the height in the site i is increased by, 

i + 1 and /ij+i > /ij + 2a and hi < /ij+2, 

i + 1 and /ij+i > hi + 2a and hi = /ii+2) 

i — 1 and > /ij + 2a and hi < /ij_2, 

i — 1 and hi-i >hi + 2a and hi = /ii_2, 



Notice that these are the rules for the growth in this column. In this context this way to describe 
the evolution equation is equivalent to derive the first moment from the master equation. In 
Fig. |l] we show this five contributions to the growth of the site i. 

The evolution equation |19| for the interface in a time step 6t = 1/N is 

hi{t + 5t) = hi{t) + 5taRi , (6) 

with 

Ri = Wi+i + W,.i + G,{h[) Wi , (7) 



1. 


a 




2. 


a/2 


if J 


3. 


a 


if J 


4. 


a/2 


if J 


5. 


a 


ifi 
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Figure 1: Schematic representation of the 5 cases that contributes to the growth of the i-th site. 
The arrows show the elected site. In grey it is shown the columns wetted at the time t. In all 
the cases the i-th. column will grow in one unit at the time t + St . 



and 

Wi±i = @{hi±i-hi-2a){[l-@{hi-hi±2)] + ShiX±2/^} , (8) 
Wi = l-@{hi-mm{hi-i,hi+i)-2a) . (9) 

Here h[ = hi + a and 6(a;) is the unit step function defined as 6(x) = 1 for x > and equals to 
otherwise. Gi{h[) equals to 1 if the cell at the height h'^ is free or active (i.e. the growth may 
occur at the next step) or if the cell is blocked or inactive. Gi is called the interface activity 
function. Taking the limit (5t ^ and averaging over the lattice we obtain {h = {hi)) 

^ = (1-W0 + (G,W0. (10) 

Here we used a = 1 as is usual in a discrete model. This equation allow us the identification of 
two separate contributions: the lateral (1 — Wi) and the local one {GiWi). 
The temporal derivative of the square interface width (DSIW) is: 

— = 2{{h-{h,))R,). (11) 

The DSIW can also be expressed by means of local and lateral additive contributions. The 
lateral contribution is 

2 [ ((1 - Wi) min(/i,_i, hi+i)) - (1 - Wi) {hi) ] , (12) 
and the local contribution is 

2[{hiGiWi)-{hi){GiWi)]. (13) 
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Figure 2: DSIW (solid line), and its lateral (Q) and local (□) contributions vs Int; for q equal 
to 0.3 (A), 0.539 (B) and 0.6 (C). 
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In order to obtain the lateral contribution [Eq. ([T2D] we used that Q{x — y) + @{y — x) — 5x,y = 1- 
In Figure |2| we plot both contributions as a function of time for various values of q. At short 
times, the lateral process is unimportant because A/i is mostly less than two. As t increases, the 
behavior of this contribution depends on q. Notice, from Eq. ([T2|), that the lateral contribution 
may be either negative or positive. The negative contribution tends to smooth out the surface. 
Figure |^ shows that this case dominate for small q. The positive lateral contribution enhances 
the roughness. This last effect is very important at the critical value. At this value, the local 
contribution is practically constant, but the lateral contribution is very strong, enhancing the 
roughness. This last contribution has important duties on the power law behavior. Generally 
speaking, the local contribution roughen the interface while the lateral one flatten it for small g, 
but the lateral contribution also roughen the interface when q increases. The lateral contribution 
is enhanced by local growth. The lateral growth may also increase the probability of local growth. 
This crossing interaction mechanism makes the lateral growth dominant near the criticality. 

6 Growth mechanisms for the Buldyrev et al. model 

The interface growth takes place under the same initial conditions as the TL model. During the 
growth, a column is selected at random with probability 1/L and the highest dry active cell, in 
the chosen column, that is nearest-neighbor to a wet cell is wetted. Afterwards, we wet all the 
dry cell below it. The height is increased by pl[ |: 

1. a if /li > max(/ij+i, and Gi{hi + a) = 1, 

2. aYi if /ij < max(/ij+i, and Gi{hi + aYi) = 1, 
where 

Y, = J2kG^{h, + ka) n {l-G,{h,+ja)) (14) 

k=l j=k+l 

and 

azi = max(/ii_i, /ij+i) - hi , (15) 

with Gi{hi + j) equal to 1 if the cell is active and if the cell is inactive. In Fig. |^ we represent 
schematically the rules. Notice that in this model the condition for the system get pinned is that 
a cluster of blocked sites, connected by nearest-neighbor expand all the lattice. This is achieved, 
in the static regime for q > qc with qc = 0.469. In this model, the time unit is defined as one 
growth attempt. In numerical simulations, at each growth attempt, the time t is increased by 
6t = 1/N. In this way, after growth attempts the time is increased in one unit. We consider 
the evolution for the height of the i-th site of the process described above. Let us denote by 
hi{t) the height of the i-th generic site at time t. Freezing the simulation at a given time, we 
compute the temporal evolution for the interface height in the next time as 

hi{t + 6t) = h^{t) + 5t a{Q{-Zi) Gi{hi + 1) + [1 - e(-z,)] Y,} , (16) 
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Figure 3: Schematic representation of the growth mechanism in the i-th site. The region in grey 
shows the height at time t, in hght grey it is shown the height at the time t + 5t. 



where Q{x) was defined in Sec.^. Yi is the increase of the height in the i-th column due to the 
contribution of the nearest-lateral-neighbor. Notice that, this kind of growth occurs by wetting 
the active cells of the chosen site nearest-neighbor of a wet cell by lateral contact, and then 
eroding all the cells bellow from a wet cell. We shall call contact contribution to the term 
[1 — Q{—Zi)] Yi and local contribution to the term B(— Zj) Gi{hi + 1) from Eq. (p^). 

Averaging over the lattice, taking 5t 0, the evolution equation for the square interface 
width is 

— = 2{{h, - {h,))Qi-z,) G,) + 2{{h, - - Qi-zi)) Yi) . (17) 

at 

The first term of both equations can be identified as the local growth contribution, and the second 
term as the contact growth contribution. In the present work we focus only on the dynamical 
behavior for the roughness. Figure |^ shows the temporal derivative of the square interface width 
(DSIW) as a function of time for various values of p. The initial condition is p in all regimes. As 



we expected [24], the power law holds only at the criticality. The DSIW goes asymptotically to 
zero at the pinning and moving phase. In Figure ^, we show the two contributions to the DSIW 
for different values of p. The local contribution 2{{hi — {hi)) 0(— Zj) Gi) to the DSIW is always 
positive. As p decreases this contribution becomes less important, but always rough the interface. 
On the other hand, for p > pc, the contact contribution 2((/ij — (hi)) (1 — Q{—Zi)) Yi) can take 
negative values, smoothing out the surface. Otherwise, for p < pc, the contact contribution 
is always positive roughening the interface. One could expect that the contact contribution 
always smooth out the surface because it tends to widen the roughen picks. However, near the 
criticality, the contact growth happens mainly in lateral neighbors cells to few height terraces 
above the mean height. Then, this new wetted column smooth out locally, but it moves away 



11 



Figure 4: DSIW as a function of time. The parameter p is 0.7 (O); 0-56 (v)) 0.531 (•) and 0.51 
(A). The symbol • shows the critical behavior. 

from the mean height increasing the roughness. 

7 Comparisons between the Buldyrev et al. and the Tang and 
Leschhorn models 

We rescue the similarities between the Buldyrev et al. and the Tang and Leschhorn models. 
In spite of the strong microscopic differences between their rules, the results obtained from 
the microscopic equation are similar. The main conclusion is that in both models the lateral 
(contact) term plays an important role in the power law behavior. The origin of the behavior 
of the lateral term near the criticality arises from a nonlinear term. This could lead to think 
that the QKPZ equation contains these features, but this is not the case as we shall show (see 
Sec. III). Moreover, the QKPZ does not allow to explain the cross mechanism between the two 
contribution, as we shall see below. 

8 Does the QKPZ describe the DPD models? 

Let us explain first why the QKPZ equation does not take into account this coupled effect. In 
the QKPZ we can distinguish two contributions, the local growth S = J-+'r](x, h) and the lateral 
one C = V d^h + ^ (dxh)"^. So, we can write the evolution equation for the height h = h{x, t) as 

il 

dth = S + C. (18) 
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Figure 5: Semi-ln plots of the different contributions to the DSIW as a function of time for the 
Buldyrev et al. model for different values of p. The circles represent the contact contribution 
and the squares represent the local contribution for the moving, the pinning and the critical 
regime, respectively. 
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Figure 6: DSIW as a function of time in the critical, pinning, and moving phases for A = 1. The 
parameter T is 0.464 (sohd line), 0.43 (dashed line), and 0.54 (dotted line). 



Taking the derivative of the square interface width, 'uP' = {{h — (/i))^), its evolution equation is 
given by 

dtw^ = 2{{h - {h))dth) = 2{{h - {h))S) + 2{{h - {h))C) (19) 

The first term can be identified as the local growth contribution, and the second term as 
the lateral growth contribution. The separation into these two analytical terms allows us to 
compare the mechanisms of growth in this model with the mechanisms in the DPD models. We 
have performed the direct numerical integration of Eq. (^) in one dimension in the discretized 
version p5|, pH| 



h{x, t + At) = h{x, t) + At {h{x - 1, t) + h{x + 1, t) 

A 9 

-2h{x, t) + - {h(x + l,t) - h{x - 1, t)Y 
8 

+T + r]{x,[h{x,t)]) } 

where [. . .] denotes the integer part and i] is uniformly distributed in [— f , f ], where a = 10^/^ is 
selected. We choose At = 0.01, use the initial condition h{x,0) = and the periodic boundary 
conditions. Figure I shows the DSIW as a function of time for various values of J-'. Here we 
found that the DSIW increases continuously from zero to a maximum value, at difference of the 
DPD models where p is the initial condition in all regimes. Moreover we did no expected to 
recover the initial regime because the QKPZ equation is valid only in the hydrodynamic limit. 



Equally to the DPD models 1 19 , |29|] , the power law holds only at the criticality. The DSIW 
goes asymptotically to zero at the pinning and moving phase. In the asymptotic regime, the 
behavior of the DSIW is similar in the QKPZ equation and in the DPD models. 

In Figure 0, we show the two contributions to the DSIW for different values of J-. The local 
contribution 2{{h — {h))S) to the DSIW is always positive. As T decreases, this contribution 
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becomes less important but always rough the interface. On the other hand, the lateral contri- 
bution 2{{h — {h))C) takes negative values in all phases smoothing out the surface. This is an 
important difference with the DPD models where, for p < p^, the lateral contribution is always 
positive roughening the interface. So, the QKPZ does not represent exactly the dynamics of 
the TL model even if the exponents derived from its numerical integration are in accord with 
the ones obtained for these models. The cross mechanism between both contributions is not 
taken into account in this equation because the noise is additive. This cross mechanism is due 
to the fact that the quenched noise is coupled to the dynamic of the interface as was shown by 
Braunstein et al. ||l9| , |20| , In order to explain the origin of this coupling let us introduce the 
passage to the continuum of our microscopic equation. 



9 Stochastic Differential Equation for the TL model 

Expanding Eq. (^) to first order in 6t, the evolution for the height of this site is 

dhi a 

— = -Ri + rji , 
ot T 



(20) 



where t = N 5t is the mean lapse between successive election of any site and rji is a Gaussian 
"thermal" noise with zero mean and covariance = (^^Z''") ^ij ^(^ ~ For 

this model. 



Ri{hi^i,hi, /ij+i) = Wi+i + Wi-i + Gi{hi + a)Wi , 



(21) 



and Wi+i, Wi-i and Wi where defined in Eq. (^). 

In order to obtain a stochastic equation we need an analytic representation of Ri. To do 
this we proceed to regularize the height defining an interpolating function |jl^ for the difference 
of height and then to expand the step function to first order in the argument as Q{x) x cq + 
cix + 0{x^). This can be done providing that x is smooth obtaining 



y ' = - [Wixi + a) + Wixi - a) + W{xi) F{x^, h{xi) + a)] + ii{xi,t) , 



(22) 



with 



W{x + a) + W{x - a) = {cq - 2ci) + A cl{d^hf + a ci 



-+4(co-2ci) 



W{x) = 1 - (co - 2ci) - 4c? {d^hf + -acidlh. 



(23) 
(24) 



Notice that the argument of G = 0(p — g{xi,h{xi) + a)) is not smooth, so its expansion is 
meaningless. In order to recover the early time regime where the dynamics is mainly random 
deposition with probability p |25, we must impose the condition cq = 2ci. The final step is 
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Figure 7: Semi- In plots of the different contributions to the DSIW as a function of time for 
different values of T and A = 1. The circles (Q) represent the local contribution, the triangles 
(A) represent the lateral contribution, and the squares (□) represent the total DSIW. The (a) 
plot shows the critical phase !F = 0.464. The (b) plot shows the pinning phase T = 0.43. The 
(c) plot shows the moving phase = 0.54. 
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a coarse-grained spatial average of the variables in order to obtain smooth continuous functions 
at a macroscopic level. In this way, we obtain the stochastic continuous equation for this model, 



where (G = G{x, h)) 



Bh 

— = ^(G) + v{G) dlh + A(G) {d^hf + 7?(x, t) , (25) 



KG) = G-, 
r 

u{G) = lci(l + G)^ 



A(G) = 4c?(l-G) 



Equation (25) shows that the nonlinearity arises naturally as a consequence of the micro- 
scopic model. As we approach to the critical value, in the dynamical regime, the density of active 
sites / = (G) goes asymptotically to zero [^] and the coefficient of the nonlinear term becomes 
relevant. In these case the main responsible of the nonlinearities is the lateral contribution [see 
Eq. ([2^)]. This explain why this contribution enhances the roughness at the criticality as was 
predicted by Braunstein et al. ||2^. Faraway above the criticality the nonlinear term becomes 
less relevant. In the limit p ^ {p ^ 1) we recover the KPZ (EW) equation with thermal noise 
as was expected. In the asymptotic regime {t ^ t*) the various derivatives of the height become 
very small on a coarse grained scale. In this temporal regime the mean height speed (MHS) 
(dh/dt) f. For p < Pc, f ^ and the MHS goes to zero while for p ^ Pc, f ^ const and 
MHS goes to constant. 

Notice that our equation is invariant under local tilting of the interface by an infinitesimal 
angle in the same way that the QKPZ equation (^) is Q, so the previous numerical results 
obtained by Amaral et al. ||T^ , that studied the effects of the effective coefficient X^jf from a 
tilted interface, are compatible with our equation. Our result is also in agreement with those of 



Reka et al. [16| that obtained numerically a parabolic shape of the local velocity as a function of 
the gradient for the DPD model near above the criticality for different reduced forces {p/pc — 1)- 
In the experiments, the advancement of the interface is determinated by the coupled effect 
of the random distribution of the capillary sizes, the surface tension and the local properties 
of the flow, so it is not surprising that all these effect give rise to a multiplicative noise. This 
multiplicative noise must be taken into account at the time to pose a model with the essential 
features of the experiment of surface growth in disordered media. In the TL and the Buldyrev 
et al. models the growing rules for the evolution of the local height are strongly coupled to the 
quenched noise in a multiplicative way. In both models the microscopic rules that allows the 



growth from an unblocked cell [19, ^] depends in some way on the local slope. In that sense this 
coupled effect is not taken into account in the QKPZ equation. The effect of a multiplicative 
noise has been proposed by Csahok et al. [^] by means of a phenomenological equation. They 
found a crossover between two temporal regimes with /3 = 0.65 to P = 0.26 but the value of 
a ~ 0.47 was obtained over a short range spatial scale. Indeed, the exponents are not the same 
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as of the DPD models. Moreover, processes with the same exponent may not belong to the 
same universahty class. For example, 1 + 1-dimensional lattice gas simulations of roughening of 
immiscible fluid-fluid interface |^] lead to the same exponents as the 1 + 1-dimensional KPZ |^] 
(/3 = 1/3 and a = 1/2) for surface growth, but this model is completely linear, so there is no 
obvious mathematical relationship between these two processes. 

10 Summary 

The microscopic equations for the models with quenched noise, treated in this work, allows to 
gain a more profound insight on the principal mechanism of growth. In special, the separation 
into the lateral and the local growth contributions allows to explain the great interplay between 
them. Obviously, the continuous equation that represents the local growth of these processes 
must take into account both mechanisms in addition of the symmetries allowed by the model, but 
this is not enough to reproduce the interplay between both mechanisms. The QKPZ equation 
contains a lateral and a local contribution but the quenched noise is additive. Notice that a 
relevant common feature of the two models treated in this work is that the rules for growth in 
active sites depend strongly on the slope. So, it is not surprising that this lead to a coupled 
effect of the quenched noise to the dynamics. This coupled effect is not taken into account 
in any equation with additive noise. Nevertheless, the QKPZ equation gives the same scaling 
exponents that these models, moreover it is not clear why. Our results suggest that the QKPZ 
equation does not describe properly the dynamics of the DPD models even if the exponents 
are similar. Our Langevin equation for the TL model reflects this coupled effect through its 
coefficient noise dependence. The equation obtained has the same terms than the QKPZ but 
its coefficient depends on the competition between the driving force and the quenched noise. In 
that sense, our equation is multiplicative in the noise. 

Finally, we conclude that the classification of these models in universality classes is not 
fully developed. The derivation of the continuous equation from the microscopic dynamics is a 
powerful method that allows to associate them in a non ambiguous way. 
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